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Abstract — Immersed boundary (IB) method has proved its 
efficacy in handling complex fluid -structure interaction 
problems in the field of Computational Fluid Dynamics 
(CFD). We present a two-dimensional computational model 
based on an immersed boundary finite volume method to 
study flow over a stationary circular cylinder. Lagrangian 
coordinates are used to describe the cylinder and Eulerian 
coordinates are employed to describe the fluid flow. The 
Navier-Stokes equations governing the fluid are solved on a 
staggered Cartesian grid system using the fractional -step 
based finite volume method. Numerical simulations are 
performed for four different Reynolds numbers to 
demonstrate the flow behavior. 

Keywords — Immersed boundary method, Flow over cylinder, 
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I. Introduction 

Fluid-structure interaction problems are most common in 
the field of computational fluid dynamics (CFD). In all the 
physical problems, the structure under consideration will be 
flexible or rigid. The flexible or rigid body exerts forces on the 
surrounding fluid and the fluid will exert forces through shear 
stresses and pressure differences. Hence, to develop a 
numerical model which can effectively describe such 
interaction to reveal the fluid dynamics involved needs careful 
attention. Conventional methods in CFD utilize body -fitted grid 
for analysing the flow around an arbitrary body. This 
methodology requires coordinate transformations and/or 
complex grid generation. For moving bodies, this technique 
further needs new mesh generation at each time, resulting in 
large computational effort. To overcome these difficulties, a 
new computational method has been developed, which is 
popularly known as the immersed boundary (IB) method. IB 
method uses a non-body conforming Cartesian grid. IB method 
is a Cartesian grid -based method in which the grid does not 
have to conform to the solid boundary. This method provides a 
modeling framework and numerical method for studying the 
interaction of incompressible fluid flow and immersed rigid or 
elastic structures. In this method, a momentum forcing function 
is incorporated in the Navier-Stokes equation to make the 
presence of the immersed boundary and the entire simulation 
can be carried out in a regular Cartesian grid. Hence the main 
advantage of IB method is that grid generation is much easier. 
Added advantages are less memory and CPU time. IB method 
was first developed by Pes kin [1] to simulate blood flow in the 
human heart. Later, various modifications have been made to 
this method by different researchers emerging it as an efficient 
numerical tool for handling complex flow simulations in the 
field of CFD. As already mentioned, the IB method utilizes 
momentum forcing function and a Cartesian grid Navier-Stokes 



solver. Based on how the momentum forcing is employed in the 
Navier-Stokes equation, IB method is classified into two: 
continuous forcing approach and discrete forcing approach. The 
continuous forcing approach of IB method utilizes Eulerian 
variables for the fluid region and Lagrangian variables for the 
solid region. The interaction between these two variables is 
linked by the Dirac delta function. The discrete version of IB 
method employs interpolation techniques to obtain the desired 
no-slip conditions at the boundary and the momentum forcing is 
obtained directly from the discretized Navier-Stokes equation. 
The continuous forcing approach is used for the simulation of 
flu id -flexible body interaction. For the simulation of elastic 
bodies interacting with the fluid, the Lagrangian force is the 
elastic force which can be obtained by applying Hooke's law. 
But when dealing with rigid bodies, this law is not well posed. 
Hence in such cases methods like the virtual boundary method 
proposed by Goldstein et al. [2] are employed. This method is 
based on the feedback forcing scheme which can enforce the no- 
slip boundary condition on the rigid boundary immersed in the 
fluid. The basic difference between the virtual boundary method 
and that of Lai and Peskin [3] is that, the boundary points are 
exactly prescribed in the former method, but allowed to move 
slightly from their equilibrium positions in the latter. Recently, 
Shin et al. [4] proposed a new version of the IB method which 
combines the feedback forcing scheme of virtual boundary 
method along with Peskin's regularized delta function. From the 
literature review on IB method, it is investigated that discrete 
forcing approach of IB method is suited for simulating rigid 
body interaction with the fluid whereas continuous forcing 
approach is suited for both rigid and elastic bodies. Also, for 
moving body problems, the continuous forcing approach of IB 
method is far superior because of the simplicity in computing the 
momentum forcing function. The present work is based on the 
continuous forcing approach of IB method. The two- 
dimensional flow past a stationary circular cylinder is a bench 
mark problem in fluid dynamics and has been widely 
investigated experimentally and numerically. Also, this problem 
has been used as the validation problem for most of the IB 
methods. With this perspective, a numerical model based on the 
Peskin's IB method [5] is developed to simulate the two - 
dimensional flow over a stationary cylinder. Here, the cylinder is 
modeled by discrete number of IB points. The Navier-Stokes 
equations governing the fluid flow over and surrounding the 
cylinder are solved on a staggered Cartesian grid system using 
the fractional-step based finite-volume method. We develop a 
FORTRAN code to perform numerical simulations to 
demonstrate the fluid dynamics of flow over a stationary circular 
cylinder. Numerical simulations are carried out for different 
values of Reynolds numbers to capture the flow behaviour over 
and surrounding the cylinder. The variation of mean drag 
coefficient with Reynolds number is also investigated. 
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The paper is arranged as follows: the solid model of the 
cylinder based on IB points and fluid model based on the 
momentum and continuity equations and the numerical 
procedure are explained in Section II. The numerical results 
generated by employing it are discussed in Section IE. Finally, 
the concluding remarks of the present study appear in Section 
IV. 

II. Mathematical Modeling and Numerical Procedure 

The general form of mathematical model based on the 
continuous forcing approach of IB method can be described as 
follows : 



— + u-Vu = -Vd + — V 2 u + f(x,r), 
dt Re v ' 



Vu = 



(1) 



(2) 



f(x,r) = \F(s,t)S(x-X(s,t))ds, 

(3) 

V(X(s,t))=\u(x,t)S(x-X(s,t))dx , 

(4) 

ax 



dt 



= V(X(s,t)) 



(5) 



where u is the fluid velocity, p the fluid pressure, Re the 
Reynolds number, f (x,^) the momentum forcing function or 
the Eulerian force density, t the time, x and s the Eulerian 
and Lagrangian variables, and <J(x)the Dirac delta function. 

Equations (1) and (2) are the Navier -Stokes and continuity 
equations for an incompressible viscous fluid in the 
dimensionless form An elastic structure or organism is 
represented as an immersed boundary at time t denoted by 

X(,s,f). Generally, the structure is assumed to have the same 

density as the fluid and hence the immersed boundary is 
assumed to be massless. The Lagrangian force F(5",f) on a 

particular segment s at time t is determined by the variational 
derivative of an elastic energy functional £"[ ] 

F(s,t) = - SE 



SX 



(6) 



The Lagrangian forces due to the structure are applied to the 
Navier-Stokes equation by means of Eq. (3). The presence of 
the structure is considered in the Navier-Stokes equation as the 

additional force f(x,^). Finally, the immersed boundary 

X^,^) is moved according to Eq. (5) using the local fluid 

velocity computed in Eq. (4). More details can be found in [1, 
3, 5]. 

The schematic view of the physical domain in the 
dimensionless form is shown in Fig. 1. Note that the non- 
dimensionalization is achieved by introducing the diameter of 

the cylinder as the characteristic length l c = d* and the inlet 
velocity of the fluid as the characteristic speed U . The 
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boundary conditions are imposed in such a way that the flow is 
from the left to right of the channel. The following boundary 
conditions are used: 

At the inlet: u = U,v = 

(7) 

At the outlet: — = 0, — = 

dx dx 



(8) 



du 



At the top wall: — = 0, v = 

ay 

(9) 



At the bottom wall: — = 0, v = 

dy 



(10) 



At the solid surface of cylinder: m=0,v = (11) 
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Fig. 1. Schematic diagram of the dimensionless physical domain 
showing the position of the cylinder. 

Here the stationary cylinder, the immersed boundary is 
considered to be elastic but extremely stiff. The cylinder is 
divided into discrete number of IB points. The Lagrangian force 
is computed using a feedback forcing method described in the 
following way, 



F(s,t) = S(X E (s,t)-X(s,t)) 



(12) 



where X(s,£) and X £ (s,t^j denote respectively the IB and its 

equilibrium location, and the stiffness constant 5»1. The 
above equation interprets that the IB points are connected to their 
corresponding fixed equilibrium points with a very stiff spring 

whose stiffness constant is S . So if the IB points fall away from 
the desired locations, the force on the spring will pull these IB 
points back. Thus, as the time goes on, it is expected that the IB 
points will always be close to their desired configurations. 

A fractional-step based finite-volume method is employed to 
solve the equations (1) and (2) to obtain the fluid velocity 



,n+l 



IT A (x) and the fluid pressure p n+ (x) at the next time level 
n + \. In this method a pseudo -pressure is used to correct the 
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velocity field so that the continuity equation is well satisfied at 
each computational time step. The governing equations are 
discretized using a third-order Runge-Kutta method (RK3) for 
the convection terms and a second -order Crank-Nicolson 
method for the diffusion terms. The step-by-step procedures for 
solving the governing equations are depicted below: 

^L = a k (L(u) + L(u k ))-2a k Vp t 
-y^(u*)- A ^(« M ) + f *' 



(13) 



vV = — *— (V-u) 

2a,At y ' 



(14) 



(15) 



u t+1 = u-2a.A?V^ 



/ +1 = /+^-^vV 

Re 



(16) 



where a k , y k di\\dp k are the coefficients of the third -order 
RK3 method, U the intermediate velocity, At the time step- 
size, (j) the pseudo -pressure introduced to satisfy the flow- 

incompressibility constraint, k the substep index and £(u) 

and Af(u)the linear diffusion and non-linear convection terms 
of the Navier-Stokes equation [6]. 

III. Results and Discussions 

The governing equations outlined in Section II are solved 
in a two - dimensional rectangular dimensionless domain of 

50 x30. The cylinder is placed at the centre of the domain. 
The Eulerian grid size is 281x221 in the streamwise (x) and 

transverse (y) directions, respectively. Eighty grid points are 

uniformly distributed inside and near to the cylinder in both the 
x- and y -directions and the grid is stretched away from the 

cylinder. The smallest grid size is 0.025. The cylinder is divided 
into 288 Lagrangian (IB) points. 

The numerical results pertinent to the flow over a 
stationary cylinder are obtained by developing a code in 
FORTRAN. Numerical simulations are performed for four 

different Reynolds numbers, Re=10, 20,40 and 100. The 

Reynolds number is based on the cylinder diameter and the 
inlet velocity. Figure 2 shows the streamlines of the flow when 
it reaches its final steady state for the three different Reynolds 

numbers ( Re=l 0,20,40). In all cases, a pair of stationary 

recirculating vortices is formed behind the cylinder. At a very 
low Reynolds number, flow around a circular cylinder is steady 
and symmetrical upstream and downstream 

As the Reynolds number increases, the upstream- 
downstream symmetry disappears and two attached vortices 
appear behind the cylinder. These vortices become bigger with 
increasing Reynolds number. As the Reynolds number becomes 
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higher than 40, alternating eddies are formed and convected in 
the wake. This generates the alternating separation of vortices, 
which are convected and diffused away from the cylinder, 
known as Karman vortex shedding phenomena. The 
instantaneous streamline pattern of unsteady laminar flow at 
Re=100 is presented in Fig. 3. The spanwise vorticity contours 
near the cylinder for Re=100 is shown in Fig. 4. The vortex 
shedding phenomena behind the cylinder is well captured in 
Fig.4. 

Table I represents the comparison of the mean drag 
coefficient obtained in the present study with those from the 
literature. In the present work, the drag forces are obtained by 
directly integrating all the momentum forcing values applied 
inside the body. A useful discussion on how to evaluate the drag 
force on a solid body in conjunction with the immersed boundary 
method can be found in Lai and Pes kin [3]. From Table I, it is 
clear that the present values of mean drag coefficient and those 
of previous research results are in good agreement. This proves 
the reliability of the developed numerical model. 

IV. Conclusion 

This paper present a two-dimensional computational model to 
study the fluid dynamics of flow over a stationary circular cylinder 
based on an immersed boundary finite volume method. 
Accordingly, the cylinder is modeled by discrete number of 
immersed boundary points. The Lagrangian forces are computed 
using a feedback forcing method. The fluid dynamics of the fluid 
flow over and surrounding the cylinder is modeled using the 
momentum and continuity equations and the set of governing 
equations are solved on a staggered Cartesian grid using the 
fractional-step based finite-volume method. 

Numerical simulations are performed using the developed 
numerical model to explore the flow behavior over the cylinder 
for four different Reynolds numbers. It is found that at low 
Reynolds numbers, flow around a circular cylinder is steady and 
symmetrical upstream and downstream For high Reynolds 
number, alternate separation of vortices, which are convected 
and diffused away from the cylinder are observed. The 
instantaneous streamline pattern for four different Reynolds 
numbers is well captured. Further, mean drag coefficients are 
computed for each case and the values are compared with the 
previous researcher's results. It is interesting to note that a good 
agreement between the present numerical results and the 
previous researcher's results is achieved. This proves the validity 
of the developed numerical model. 
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Fig. 2. Streamline visualization of the steady flow over a circular cylinder at 
Re=10, 20 and 40 
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Fig. 4. Vorticity contours nearthe circular cylinder for Re=100 

TABLE I 

Comparison of mean drag coefficient for different Reynolds numbers 



Reynolds 


Su et al. 


Silva et al. 


Deng et al. 


Present 


Number 


[7] 


[8] 


[9] 




10 




2.81 


2.98 


2.98 


20 


2.20 


2.04 


2.06 


2.16 


40 


1.63 


1.54 


1.52 


1.67 


100 


1.40 


1.39 


1.30 


1.32 



Fig. 3. Instantaneous streamline pattern of the unsteady flow observed at 
Re=100 
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